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A  fully  Sinc-Galerkin  method  for  recovering  the  spatially  varying  stiffness  parameter 
in  fourth-order  time'dependent  problems  with  fixed  and  cantilever  boundary  conditions  is 
presented.  The  forward  problems  are  discretized  with  a  sine  basis  in  both  the  spatial  and 
temporal  domains.  This  yields  an  approximate  solution  which  converges  exponentially  and 
is  valid  on  the  infinite  time  interval.  When  the  forward  methods  are  applied  to  parameter 
recovery  problems,  the  resulting  inverse  problems  are  ill-posed.  Tikhonov  regularization  is 
applied  and  the  resulting  minimization  problems  are  solved  via  a  quasi-Newton/trust  region 
algorithm.  The  L-curve  method  is  used  to  determine  an  appropriate  value  of  the  regulariza¬ 
tion  parameter.  Numerical  results  which  highlight  the  method  are  given  for  problems  with 
both  fixed  and  cantilever  boundary  conditions.  ^ 
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1  Introduction 


In  this  paper,  a  fully  Sinc-Galerkin  method  is  introduced  for  the  numerical  recovery  of  mate¬ 
rial  parameters  in  fourth-order  time-dependent  problems.  To  illustrate  the  method,  consider 
the  problem  of  estimating  the  spatially  varying  parameter  EI(x)  in  the  state  equations 


<  x  <  1  t  >  0 


u(0,£)  =  tt(l,f)  =  0,  £>0 

|(0,<)  =  |(M)  =  0,  *>0 

Q 

u(x,0)  =  -^(x,0)  =  0)  0  <  x  <  1 


(1.1) 


and 


C(EI)u  =  /(x,  £),  0  <  x  <  1,  t  >  0 

u(0,t)  =  a(0.  (b/^)(1,0=T(<).  i>0  (12) 

=  |(£/g)(M)  =  J(0,  ‘>o 

U(X,  0)  =  =0>  0  <x  <  1 

given  measurements  of  the  data  at  the  points  {(xp,  tq  )}p=l[  in  (0, 1)  x  2R+.  These  formula¬ 
tions  are  generalizations  of  the  equations  which  arise  when  using  the  Euler-Bernoulli  theory 
to  model  beams  with  flexural  rigidity  EI(x)  and  fixed  and  cantilever  ends,  respectively.  For 
ease  of  presentation  throughout  the  paper,  the  boundary  conditions  in  (1.1)  and  (1.2)  will 
be  referred  to  as  fixed  and  cantilever  conditions  with  the  general  7(f)  and  £(f)  included  to 
allow  for  boundary  controllers. 

Since  EI(x)  denotes  the  flexural  stiffness,  it  is  physically  reasonable  to  assume  that  El 
is  continuous  on  [0, 1]  and  to  let  the  admissible  parameter  set  Q  be  defined  by 


Q  =  {EIe  H\ 0, 1) :  EI(x)  >  EI0  >  0} 


(see  (9)).  With  this  definition,  the  existence  of  a  unique  solution  u  to  the  forward  problem 
can  be  obtained  on  any  fixed  time  interval  [0,  r],r  >  0,  for  /  sufficiently  smooth. 
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In  order  to  apply  the  results  from  classical  operator  theory,  the  inverse  problems  corre¬ 
sponding  to  (1.1)  and  (1.2)  can  be  written  as  operator  equations  of  the  form  K(EI)  =  d 
where  K  is  compact.  The  procedure  for  problems  with  fixed  boundary  conditions  is  out¬ 
lined  below  with  the  formulation  for  problems  with  cantilever  boundary  conditions  following 
similarly.  Further  theory  for  this  latter  case  can  be  found  in  [1]. 

In  formulating  the  operator  equation  corresponding  to  (1.1),  it  is  easiest  to  first  consider 
the  spatial  problem 


(EI(x) u")"  —  q(x)u  =  /(x),  0  <  x  <  1 

u(0)  =  u(l)  =  u'(0)  =  u'(l)  =  0 


(1.3) 


where  EI{x)  and  g(s)  are  both  strictly  positive  on  [0,1].  From  the  theory  of  [10]  as  noted 
in  [20],  there  exists  a  Green’s  function  for  (1.3)  which  will  be  denoted  by  G{xtstEI(a))  to 
emphasize  its  dependence  on  El.  It  follows  that  G  and  are  continuous  on  [0,1]  and  that 
the  state  solution  u  is  given  by 


tt(x)  =  f  G(x,s,EI(s))f(s)ds. 
Jo 


(see  [20],  pages  205-206). 

These  one- dimensional  results  can  then  be  extended  to  the  time-dependent  problem  (1.1) 
via  a  separation  of  variables.  The  substitution  of  u(x,t)  =  X(x)T(t)  into  the  homogeneous 
problem  corresponding  to  (1.1)  yields 


and  the  eigenvalue  problem 


T'\t)  +  XT{t)  =  0 


(EI(x)X"(x))"  -  \X(x)  =  0,  0  <  *  <  1 

*(0)  =  JT(1)  =  *'(0)  =  .Y'(l)  =  0 


where  A  >  0.  From  arguments  similar  to  those  in  [5]  and  [20],  it  follows  that  since  EI(x)  >  0 
on  [0,1],  the  state  solution  to  (1.1)  can  be  represented  by 


The  function  Q  depends  on  a  Green’s  function  as  well  as  the  expansions  of  T(t)  and  the  forcing 
function  /.  Hence  Q  is  continuous  with  respect  to  El  and  is  bounded  on  (0, 1)  x  (0,r] 
for  /  sufficiently  smooth. 

The  inverse  problem  can  then  be  written  as  the  operator  equation 

IC(EI)  =  d  (1.4) 

where  d  denotes  the  data.  The  nonlinear  operator  K, :  H2( 0, 1)  — *  La((0, 1)  x  (0,  r])  is  defined 

by  ^ 

K.(EI)  =  C  [l  G(’,-ts,EI(3))d3  (1.5) 

Jo 

where  the  observation  operator  C  maps  the  state  solution  into  the  infinite  dimensional  ob¬ 
servation  space  by 

=  {#M)>  ;  (1-6) 

that  is,  C  samples  functions  continuously  throughout  (0,1)  x  (0,r).  Note  that  the  choice 
LJ((0, 1)  x  (0,  r])  for  observation  space  is  physically  reasonable. 

To  show  that  K  is  compact,  it  is  first  shown  that  it  is  weakly  continuous;  that  is, 

EIn  El  in  H2( 0, 1)  implies  that  K(EIn )  -♦  K(EI)  in  the  L 2  norm.  For  every  (x,  t) 

in  (0, 1)  x  (0,  r]  it  follows  that 

I cf  G(x,t,3,EIn(3))ds-C  I  G(x,t,s,EI(s))ds 
I  Jo  Jo 

<  l'  \EIn(s)  -  EI(,)\i, 

<M  I'  |E/„(j)  -  EI{s)\ds 
Jo 

<  M  max  |£/„(3)  -  EI(s)\ 

with  the  last  inequality  resulting  from  the  continuity  of  El  on  [0,1].  Since  weak  conver¬ 
gence  in  J/a( 0, 1)  implies  uniform  convergence,  it  follows  that  fC(EI„)  converges  uniformly 
to  JC(EI).  The  weak  continuity  of  K  results  from  the  fact  that  the  uniform  norm  is  stronger 
than  the  1 2  norm.  Hence  the  operator  K  as  defined  in  (1.5)  is  compact  since  weak  continuity 
implies  compactness. 

Although  the  procedure  just  outlined  was  for  the  problem  (1.1),  similar  results  can  be 
obtained  for  (1.2)  once  a  Green’s  function  has  been  found  which  satisfies  the  boundary 
conditions.  Further  theory  on  problems  of  this  type  can  be  found  in  [1]  and  [10]. 


3 


Consider  now  the  well-posedness  of  (1.4).  First,  since  Q  is  continuous,  the  range  of 
the  operator  K  lies  in  C7((0, 1)  x  (0,r])  for  any  EI  €  Q.  Hence  there  exist  elements 
d  €  L2((0, 1)  x  (0,r])  for  which  (1.4)  has  no  solution.  Furthermore,  since  K  is  a  compact 
operator  with  an  infinite  dimensional  range,  it  follows  that  the  Moore-Penrose  generalized 
inverse  £*  is  discontinuous.  This  in  turn  indicates  that  small  perturbations  in  the  data  d 
may  give  rise  to  arbitrarily  large  perturbations  in  the  solution  El  €  Q.  Consequently,  some 
sort  of  regularization  (i.e.,  stabilization)  is  required  to  obtain  an  accurate  approximation  for 
EI. 

The  regularization  technique  that  is  used  is  Tikhonov  regularization  [24] ,  and  the  problem 
(1.4)  is  replaced  by  the  minimization  problem 

minr„(£7)  (1.7) 

where 

ubi)  =  i{n  mi)  -  41’ + 

Here  a  >  0  is  a  regularization  parameter  which  controls  the  tradeoff  between  goodness  of  fit 
to  the  data  and  stability.  The  penalty  functional  J{EI)  provides  stability  and  all  ws  the 
inclusion  of  a  priori  information  about  the  true  parameter  EI.  Since  EI  is  assumed  to  be 
“smooth”  in  the  sense  that  EI  €  f/2(0, 1),  the  penalty  functional  is  taken  to  be  the  norm 

J(E1)  =  ||£/||»  =  jf'[£/"(x)]Vi  +  c^'[E/(i)]Vx.  (1:8) 

with  e  of  order  10-a.  The  reasons  for  including  the  second  term  and  forcing  J  to  be  strictly 
positive  will  be  discussed  in  the  fourth  section  of  the  paper.  By  using  arguments  similar  to 
those  in  [7]  and  [16]  and  assuming  that  IC(EI)  is  one  to  one,  it  can  be  shown  that  with  this 
definition  for  J(EI),  the  solutions  EIa  to  (1.7)  converge  as  the  regularization  parameter 
a  — »  0  and  as  the  perturbations  in  the  data  and  operator  tend  to  zero. 

Due  to  the  infinite  dimensionality  of  Q  and  that  of  the  state  space,  the  problem  (1.7)  is  an 
infinite  dimensional  minimization  problem.  In  order  to  develop  a  practical  numerical  scheme, 
the  problem  must  be  replaced  by  a  sequence  of  finite  dimensional  problems;  that  is,  one 
must  approximate  the  operator  K  and  minimize  the  functional  Ta  over  a  finite  dimensional 
admissible  subspace  of  Q. 
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The  evaluation  of  tC(EI)  requires  the  solution  of  the  partial  differential  equations  (PDE’s) 
(1.1)  or  (1.2).  Similar  PDE’s  must  be  solved  to  obtain  the  components  of  the  derivative 
tC'(EI).  The  construction  of  an  approximate  solution  to  these  forward  problems  commonly 
begins  with  a  Galerkin  discretization  of  the  spatial  variable  with  time-dependent  coefficients. 
This  yields  a  system  of  ordinary  differential  equations  which  is  solved  via  differencing  tech¬ 
niques.  Due  to  stability  constraints  on  the  discrete  evolution  operator,  low  order  methods 
with  small  time  steps  are  often  required  to  obtain  accurate  approximations.  Moreover,  this 
time-stepping  must  be  repeated  at  each  step  in  the  minimization  of  (1.7).  A  final  difficulty 
lies  in  the  need  to  interpolate  at  data  points  which  do  not  coincide  with  the  nodes  of  the 
ODE  solver. 

In  contrast,  the  method  of  this  work  implements  a  Galerkin  scheme  in  time  as  well  as 
space.  This  method  thus  bypasses  many  of  the  difficulties  associated  with  time-stepping 
methods  in  the  context  of  inverse  problems.  Corresponding  results  for  the  heat  equation  can 
be  found  in  [12]  and  [19]. 

The  fully  Sinc-Galerkin  method  in  space  and  time  has  many  salient  features  due  both 
to  the  properties  of  the  basis  functions  and  the  manner  in  which  the  problem  is  discretized. 
Perhaps  the  most  distinctive  feature  of  the  method  is  the  resulting  exponential  convergence 
rate  when  solving  the  corresponding  forward  problems.  Furthermore,  the  judicious  choice 
of  a  conformal  map  provides  approximate  solutions  which  are  va’id  on  the  infinite  time 
interval  rather  than  only  on  a  truncated  time  domain.  Finally,  the  discrete  system  requires 
no  numerical  integrations  to  fill  either  the  coefficient  matrices  or  the  right-hand  side  matrix. 
All  three  features  prove  to  b*;  advantageous  when  solving  the  forward  problems  and  hence 
the  inverse  problem. 

The  foundations  of  the  Sinc-Galerkin  method  are  described  in  Section  2.  The  fundamental 
quadrature  rules  are  given,  and  the  exponential  convergence  rate  of  this  method  is  stated. 
A  thorough  review  of  sine  function  properties  can  be  found  in  [22]  and  [23]. 

In  the  third  section,  the  Sinc-Galerkin  systems  for  the  forward  problems  are  constructed 
and  implementation  details  are  discussed.  The  section  includes  the  outline  of  a  very  robust 
and  accurate  algorithm  for  solving  the  resulting  matrix  systems. 

Section  4  includes  the  finite  dimensional  minimization  problem  with  the  discussion  cen- 
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tering  around  the  construction  of  the  various  components  of  the  Tikhonov  functional.  The 
resulting  unconstrained  optimization  problems  are  solved  via  a  quasi-Newton/trust  region 
algorithm  as  described  in  [2]  and  [8]. 

Numerical  results  are  presented  in  Section  5.  Of  the  many  examples  tested,  those  dis¬ 
cussed  in  this  section  best  exhibit  the  features  necessary  for  the  practical  implementation 
of  the  method.  A  brief  discussion  of  the  Z-curve  technique  [6]  for  determining  the  regular¬ 
ization  parameter  a.  is  given  at  the  beginning  of  the  section,  and  the  applicability  of  this 
technique  in  conjunction  with  the  Sinc-Galerkin  method  is  demonstrated  by  the  numerical 
results.  Finally,  results  are  included  both  from  data  sets  with  white  noise  and  from  data  sets 
to  which  no  noise  was  added.  As  shown  in  these  examples,  the  Sinc-Galerkin  method  works 
equally  well  in  both  cases. 


2  Sine  Function  Properties 


For  the  Sinc-Galerkin  method,  the  basis  functions  are  derived  from  the  Whittaker  cardinal 
(sine)  function 


stnc(x) 


8in(irx) 

— - - ,  — oo  <  x  <  oo 

irx 


(2.1) 


and  it  translates 


5(fc,  h)(x)  =  sine 

For  h*  =  three  adjacent  members  of  this  sine  family 
in  Figure  1. 


,  h  >  0. 

(S(fc,h*)(x),fc  =  —1,0,1)  are  shown 
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x-xxU 

Figure  1.  Three  Adjacent  Members  (S(k,  h*)(x),k  =  —1,0,1,  h*  =  of  the  Translated  Sine 
Family. 

To  construct  basis  functions  on  the  intervals  (0,1)  and  (0,oo),  respectively  consider  the 
conformal  maps 

«*)  =  In  (j^)  (2.2) 

and 

T(u>)  =  ln(w).  (2.3) 

The  map  <f>  carries  the  eye- shaped  region 

DB  =  {z^x  +  iy:\arg[T^)\<d<\)  (2.4) 

onto  the  infinite  strip 

0s  =  K  =  <  +  "rlil<<*<f}-  (2.5) 

Similarly,  the  map  T  carries  the  infinite  wedge 

Dw  =  {u>  =  t  +  is  :  |arp(u/)|  <  d  <  ~}  (2.6) 

onto  the  strip  Ds .  These  regions  are  depicted  in  Figure  2. 
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Figure  2.  The  Domains  Ds,  De,  and  Dw- 


The  sine  gridpoints  Zk  €  (0, 1)  in  De  will  be  denoted  x*  since  they  are  real.  Similarly,  the 
gridpoints  Wk  €  (0,  oo)  in  Dw  will  be  denoted  t*.  Both  are  inverse  images  of  the  equispaced 
grid  in  Ds',  that  is, 

ekh 

xk  =  r\kh)  =  (2.7) 

and 


tk  =  r~\kh)  =  ekh. 


(2.8) 


To  simplify  notation  throughout  the  remainder  of  this  section,  the  pairs  De  and  T,Dw 
are  referred  to  generically  as  x»  D.  It  is  understood  that  the  subsequent  definition,  theorems, 
and  identities  hold  in  either  setting.  Furthermore,  the  inverse  of  x  is  denoted  by  ip. 

The  important  class  of  functions  for  sine  interpolation  and  quadrature  is  denoted  B{D) 
and  defined  next. 
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Definition  2.1.  Let  B{D )  be  the  class  of  Junctions  F  which  are  analytic  in  D,  satisfy 


f  \F(z)dz\  -+0,  t  — ►  ±oo 


where  L  =  {ts  :  |s|  <  d  <  and  on  the  boundary  of  D  (denoted  dD)  satisfy 


ff(F)  5  /  \F(z)dz\  <  oo. 
JBD 


The  following  theorem  for  functions  in  B(D)  is  found  in  [21]. 

Theorem  2.1.  Let  T  be  (0,1)  or  (0,oo)  when  \  =  <f>  or  T,  respectively.  If  F  €  B(D)  and 
Zj  =  ip{jh)  =  x~lWt  j  =  0,  ±1,  ±2,  •  •  • ,  then  for  h  >  0  sufficiently  small 


X  n*)k  -  ^  ^ 


Theorem  2.1  illustrates  the  exponential  convergence  rate  which  is  a  trademark  of  the  sine 
methods.  There  is  a  common  occasion  when  it  is  possible  to  evaluate  the  infinite  series 
appearing  in  (2.9),  namely  when  integrating  against  S(k,h )  o  In  general,  however,  the 
series  must  be  truncated.  With  additional  hypotheses,  it  is  proven  in  [11]  and  [22]  that  the 
truncation  need  not  be  at  the  expense  of  the  exponential  convergence. 

Theorem  2.2.  Assume  F  €  B(D)  and  that  there  exist  positive  constants  K,a,  and  f3  such 


F(t )  I  «  °,X(T)I,  rG^((-oo,0)) 

x'(T)  ~  K  \ 

I  e“^WTH,  r  €  ^([0,00)). 


(2.10) 


Then  for  h  sufficiently  small 


X  + z + f . 


(2.11) 


Theorems  2.1  and  2.2  are  used  to  establish  a  uniform  error  bound  when  constructing  an 
approximate  solution  to  the  forward  fourth-order  time-dependent  problems.  The  application 
of  these  quadrature  theorems  is  facilitated  by  the  identities 


=  (•S'Cp,  0  x(«)]| 

l*=* 


1,  i  =  p 

0,  i^P, 


(2.12) 
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4*1}  =  h  4-S(P>h)ox(z) 

•  ^  *=*4 

-Tl S(p}h)oX{z)  =| 

-  ^  f=fi  I 


0, 

(-iy-p 

(*-p) 


* =p 


»  »7^P. 


3,  *  =  P 

(_2)(_1)*-p 

i  =  p 

•[6  —  ira(i  —  p)a]»  i/p, 


4-’  =  h’  5(P> A) 0  x(z)  ^ 

I  (•  -P)1 


(2.13) 


(2.14) 


(2.15) 


=  l*'  ^S(p,h)oX(z) 


|  -fr-p/l6  "  »’(<  -  l’)1)-  '  *  P 


(2.16) 


which  denote  the  evaluation  at  the  gridpoint  Z;  of  the  sinc-map  compositions  and  their 
derivatives  with  respect  to  the  map 


3  The  Forward  Problem 


Two  forward  problems  of  interest  are 


giu  fii  /  Q2U  \ 

Cu(x,t)  =  ^-(x.t)  +  fat  ^EI(x)^(xtt)J  =  /(*»*)»  0  <  * 


<  1  <>0 


u(0,t)  =  u(l,t)  =  0,  <>0 

£(0.0- 5(1,0 -0.  *>0 


u(x,0)  = -^(x,0)  =  0,  0  <  X  <  1 
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and 


Cu(x,t)  =  /(®,t),  0<*<1,  t>0 
u(0, ()  =  *(«).  (e/^)(1,0  =  T(0.  ‘>0  (3.2) 

£(0.o -m  s(«£) (>.«)-*( 0.  ‘>0 

du 

u(x,0)  =  -^(x,0)  =  0,  0  <  x  <  1  . 

Since  a  thorough  derivation  of  the  Sinc-Galerkin  method  for  problems  of  this  type  is 
given  in  [18],  the  following  discussion  contains  only  that  material  which  is  needed  for  the 
construction  of  the  associated  matrix  systems. 

To  define  the  Sinc-Galerkin  approximation  to  (3.1),  let  5<(x)  =  S(i,hx )  o  $(x)  and 
S?(f)  =  S(j,  ht)  o  T(t),  and  take  the  basis  to  be  {  where 

“  Si^x'jSj (t). 

The  approximate  solution  is  then  defined  by  way  of  the  tensor  product  expansion 

at.  ff,  mx  =  Mx  +  Nx  + 1 

0  =  53  5-/  UijFij(x,t),  (3.3) 

i=-M.  i=-M,  m,  =  Mt  +  Nt  +  1 . 


The  mx  •  mt  unknown  coefficients  {tt,y}  are  determined  by  orthogonalizing  the  residual  with 
respect  to  the  set  of  sine  functions  {SpS*}pZT m.’,.  ’, •  This  yields  the  discrete  Galerkin 
system 


{Cummmt  -  f,SpS*)  =  0 

(3.4) 

for  p  =  —  M*,  •  •  • 

,  Nx  and  q  =  -Mt,  •  •  • ,  Nt.  The  inner  product  (•,  •)  is  taken  to  be 

(F,G)=  /  /  F(x,t)G(x,t)w(x,t)dxdt 

Jo  Jo 

(3.5) 

with  the  weight 

w(x,t)  =  w(x)w\t)  =  (f  (i))-»(f  (())-». 

(3.6) 
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The  expressions  (3.1),  (3.4),  and  (3.5)  are  combined  to  form  the  system 

+  Jo  J fQ  ^  (EI(x)^umm>mXxtt)j  Sp{x)Sl{t)w(x)w*(t)dxdt  (3.7) 

=  f  f  f(x,t)Sp(x)S*(t)w(x)w*(t)dxdt 
Jo  Jo 

for  p  =  —Mxt  ■  •  • ,  Nx  and  q  =  —Mt,  •  •  • ,  Nt. 

In  anticipation  of  the  parameter  identification  problem  which  motivates  this  analysis,  the 
term  EI(x)  in  (3.7)  is  expanded  as  a  linear  combination  of  weighted  sine  functions  with  four 
Hermite-like  algebraic  terms  added  to  accommodate  the  potentially  nonzero  function  and 
derivative  values  of  El  at  *  =  0  and  x  =  1.  Specifically,  this  parameter  basis  is  taken  to  be 
with 

b-Mm{x) ,  k  =  —Mt 

6_m.+i(*)i  k  —  -Mx  +  1 

Mx )  =  uB(x)S*(x),  Jfc  =  -M,  +  2,  •  •  ■ ,  Nx  -  2  (3-8) 

bffm- i(*)>  k  =  Nm-1 

^jv,(a:)»  k  =  Nx. 

i 

Here  S*(x)  =  S(k,hx)  o  <f>{x )  and  the  basis  weight  function  vE  is 

vE(x)  =  w(x)  —  [x(l  —  x)]s.  (3.9) 

The  algebraic  boundary  basis  functions  are  given  by 

&-M.+ 1(*)  =  (1  ~  ®)J[2®  +  1], 

bNm- iix)  =  *2[2(1  —  *)  +  1], 

)  =  ®(1  -  *)’, 

and 

bN.(x)  =  “(I  -  x)x *• 
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(3.10) 


The  finite  dimensional  approximation  of  El  then  takes  the  form 

EIma(X)  =  12  c*V»fc(*) 

k=-Ma 

where  the  choices  Me  =  Mx  and  Ne  =  Ns  are  made  so  as  to  guarantee  a  square  spatial 
coefficient  matrix.  In  the  forward  problem,  the  coefficients  {c*}£f_Ma  are  known  whereaB  in 
the  corresponding  parameter  recovery  problem,  they  are  unknown  and  are  determined  via 
methods  to  be  discussed  in  Section  4. 

A  quick  note  should  be  made  concerning  the  choice  of  parameter  basis  and  the  manner 
of  expanding  EIma.  The  two  derivative- interpolating  boundary  basis  functions  are  added 
so  that  this  expansion  of  EIma  is  the  same  as  that  used  with  cantilever  or  free  boundary 
conditions.  The  choice  of  (3.9)  for  basis  weight  is  certainly  sufficient  and  proves  to  be 
beneficial  when  incorporating  this  forward  scheme  into  a  numerical  method  for  solving  the 
parameter  recovery  problem  as  described  in  Section  4. 

The  expansion  (3.10)  is  substituted  into  (3.7),  and  integration  by  parts  is  used  to  transfer 
the  derivatives  onto  the  product  SpwS*w*.  As  detailed  in  [18],  the  weight  choice  (3.6)  guar¬ 
antees  that  all  boundary  terms  vanish.  The  resulting  integrals  are  evaluated  via  Theorem  2.2, 
or  when  possible,  Theorem  2.1.  The  requirement 

|£J(x)u(x,  t)|  <  /fxa+*(l  —  x)^+^<'T+^f*  , 
where  the  “homogeneous”  part  of  El  is 

61{x)  =  EJ(x)  -  EI(0)b_Mm+1(x)  -  £/(l)&*._i(*)  -  Er(0)6_Af.(x)  -  £J'(l)&*.(x),  (3.11) 

guarantees  the  decay  needed  to  truncate  the  infinite  quadrature  rule  as  specified  by  (2.10). 
With  a,/?, 7  and  S  specified  and  Mx  chosen,  the  choices 


/  ird 

(3.12) 

hm  - 

h 

t  =  hmi 

(3.13) 

N.= 

[>+1]- 

(3.14) 

M,  = 

2m. +  1]. 

(3.15) 
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and 


N,  =  []M,  +  l]  (3.16) 

for  the  stepsizes  and  summation  limits  balance  the  asymptotic  errors  to  at  least  order 
0  .  This  rate  results  from  the  presence  of  a  sine  function  in  the  integral.  In  the 

above  expressions,  [•]  denotes  the  greatest  integer  function.  Note  that  the  +1  is  unnecessary 
when  |Af„  ^Mxt  or  ^Mt  are  integers. 

In  many  time-dependent  problems,  the  solution  decays  exponentially  at  infinity;  that  is, 
the  solution  satisfies 

\El(x)u{x,t)\  <  tfxa+5(  1  -  xyMV+*e-w.  (3.17) 

With  this  supposition,  Lund  [11]  shows  that  the  condition  (3.16)  can  be  replaced  by 

*  =  lh,n  (?MA) + ■]  •  <318> 

The  selection  Nt  in  (3.18)  significantly  reduces  the  size  of  the  discrete  system  with  no  loss 
of  accuracy. 

Given  Ms,  NSt  Mt,  Nt  and  h  =  hm  =  ht  as  defined  above,  the  discrete  system  for  (3.1)  is 
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It  should  be  noted  that  the  ordering  of  the  coefficients  tt,j  in  U  mimics  that  used  in  most 
standard  time-differencing  schemes.  This  is  a  matter  of  convenience  since  the  Sinc-Galerkin 
method  is  not  bound  by  any  specific  ordering  of  the  grid. 

To  simplity  notation  when  specifying  the  spatial  and  temporal  matrices  A(EI)  and  At 
respectively,  let  /W,  l  =  0,1,2, 3,4  denote  the  matrices  whose  pt-th  entry  is  from 
(2.12)  -  (2.16).  As  shown  in  [13],  the  mt  x  mt  matrix  At  is  given  by 

At=  U((f)*).  (3.23) 

The  m«  x  mx  matrix  A(EI)  has  the  form 

A(E1)  =  [*<J)2?(p*<„)  +  2*<3>2>(j?*(0)  +  *(4>2>(p*<o>)] .  (3.24) 

The  notation  P(p*(o)>  t  —  0, 1,2  denotes  the  diagonal  matrices  containing  the  components 

of  the  vectors 

pm  =  t  =  0,1,2  (3.25) 

where  c—  [c_m., •  •  •  ,cjv.]T-  The  matrices  j  =  2,3,4  and  l  =  0,1,2  are  defined 
componentwise  by 

I*0!'  =  ^(VPM  (3.26) 

and 

(*«)),. t  (3.27) 

with  the  notation  on  the  right-hand  sides  of  (3.26)  and  (3.27)  indicating  the  j- th  and  /-th 
derivatives,  respectively. 

To  illustrate  the  dependence  of  4>W,  j  =  2,3,4  and  l  =  0,1,2  on  fundamental 
matrices,  the  respective  expansions  are  listed  below.  The  diagonal  matrices  V  and  the 
matrices  /W,  l  =  0,1,2, 3, 4  have  sizes  consistent  with  the  following  range  of  indices  »,p, 
and  k  (—A/*  <  t  <  NXi  —Mx  <  p  <  NXi  —Mx  -f  2  <  k  <  Nx  —  2).  From  (3.26)  it  follows  that 
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and 


$<»)  =  JL/(»P(ti,^)  +  +  7<°>X>(  w"), 

hi  ** 


$(3)  =  _L/(3)£>  (u;(^')a)  +  ~J(2)I?(3tyy  +  3«)f) 

+EJW®  (3“" + 3wj + •?)  +  ,Wt>  (?) 


+ 

+ 


^l{7)v  [tor 4/  +  12t u'f  +  4 w<f>'"  + 


j-/(1)p  ( 4ii/"  +  6u>"^  +  4u>'^  +  ti;— ) 
hx  \  <p'  <(>  <f>  J 

(t) 


+/(°)p 


For  l  =  0, 1,2  the  mx  x  mx  matrices  in  (3.27)  are  given  by 

*«  =  ;  fl(,) :  JSL, :  # 


(3.28) 

(3.29) 

(3.30) 


(3.31) 


where  6^  =  •  •  • , t>k\xN.)]T  for  fc  =  -Mx, -  AT*  +  1,  JVX  -  1  and  Nx.  Again,  the 

superscript  l  indicates  the  £-th  derivative.  The  mx  x  (mx  —  4)  matrices  are 


BW  =  V(vB)lW,  (3.32) 

BW  =  'v(veV)IM  +  V(vE)lW  (3.33) 

and 

B(2)  =  tfV  {M<t>')3)  *{a)  -  Yv(VE<t>' +  2u^')/(1>  +  p(^)/(0)-  (3-34) 

The  negative  signs  that  appear  in  the  definitions  of  B W  and  B ^  result  from  the  transposing 
of  Again,  it  is  noted  that  in  (3.32)  -  (3.34),  the  mB  x  (m,  —  4)  matrices  I^l\t  =  0, 1,2 
have  components  as  defined  in  (2.12)  -  (2.14). 

Various  methods  exist  for  solving  matrix  systems  of  the  form  (3.19),  one  of  which  derives 
from  the  generalized  Schur  decomposition  (page  396  of  [4]).  As  guaranteed  by  the  results  of 
Moler  and  Stewart  [15],  there  exist  unitary  matrices  Qi,ZltQ3  and  Z3  such  that 
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Q\AXZX  =  P 
Q\CmZ\  =  R 


Q*2CtZ2  —  S 


Q\AtZ2  =  T 

where  P,R,S,  and  T  are  upper  triangular.  If  Y  =  Z[UZ2  and  C  =  Q\GQ2 ,  then  (3.19) 
transforms  to 


PYT*  +  RYS*  =  C. 


By  comparing  the  /b-th  columns,  one  finds  that 

n  n 

P  tkiVi  +  R  Y,  skiVj  =  cfc 

j=k  j=k 

which  yields 

n  n 

(tkkP  4-  &kkR)yk  =  c*  —  P  JZ  tkiVi  —  R  akjy>  (3.35) 

ji=fc+i  >=fc+ 1 

(for  convenience,  it  is  assumed  that  all  matrices  are  n  x  n  and  indexed  from  1  to  n).  With 
the  assumption  that  the  matrix  ( +  SkkR)  is  nonsingular,  the  solution  to  (3.35)  is  easily 
found  by  recursively  solving  triangular  systems. 

Although  this  algorithm  does  require  complex  algebra,  it  is  both  robust  and  efficient 
and  requires  no  assumptions  concerning  the  diagonalizability  of  the  component  matrices. 
It  should  be  noted  that  a  “real”  version  of  this  algorithm  also  exists  [3].  In  this  latter 
algorithm,  Qi,Zi,Q2,  and  Z2  are  orthogonal  with  P,S  quasi-upper  triangular  and  R,T 
upper  triangular. 

A  Sinc-Galerkin  method  for  the  more  general  problem  (3.2)  can  be  developed  in  a  similar 
manner  once  a  suitable  basis  has  been  determined  for  discretizing  the  spatial  variable.  To 
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this  end,  define  the  set  of  spatial  basis  function  by 


w  = 


B-M.- 4(x), 

i  =  -Ms  -  4 

B-m.- a(*), 

•*. 

II 

1 

.* 

1 

CO 

B-m.- 3(x), 

i  =  -M.  -  2 

B-m.- i(x), 

t  =  -Mm  -  1 

v(x)Si(x), 

t  =  —  Afe,  •  •  • 

Bn.+ i(*)» 

*  =  N,  +  1 

bn.+ a(x), 

t  =  JV.  +  2 

Btf.+ s(x), 

*  =  ^,  +  3 

Bn.+4(*), 

i  =  Nx  +  4. 

Here  Si{x )  =  S(i,  hx)  o  <f>[ x)  and  the  basis  weight  v(x)  is  taken  to  be 

v(x)  =  [x(l  -  x)]3. 

The  boundary  basis  functions  are 

=  (1  -  x)4[20x3  +  10x2  +  4x  +  1], 

BNm+i(x)  =  x4[20(l  -  x)3  +  10(1  -  x)2  +  4(1  ~  x)  +  1], 
B-m.- a(x)  =  x(l  —  x)4[10x2  +  4x  4-  1], 

Bn.+ a(x)  =  x4(l  -  x)[10(l  -  x)2  +  4(1  -  x)  +  1], 
B.u..,(x)  =  x’(l  -  ,)*  [ix  +  j] , 

Bw.+i(x)  =  x*(l  -  xf  [2(1  -  *)  +  i]  , 


and 


B-m.- <(*)  =  ^*3(1  -  x)4, 


Bn.+a(x)  =  -ix4(l  -  x)3. 


(3.36) 


(3.37) 
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A  brief  note  concerning  the  choice  of  spatial  basis  is  in  order  at  this  point.  First,  since 
hm)  o  ^(x)]  ,  l  =  1, 2,  •  •  • ,  are  undefined  at  x  =  0  and  x  =  1,  some  basis  modifica¬ 
tions  must  be  made  when  solving  problems  with  nonzero  boundary  conditions  (see  also  the 
definition  of  ipk  in  (3.8)).  With  the  cantilever  boundary  conditions  of  (3.2),  it  is  tempting  to 
use  fewer  algebraic  boundary  basis  functions  and  the  basis  weight 

u(x)  =  ®(1  -  x)3 


but  in  many  problems  this  results  in  nonzero  boundary  terms  when  integrating  by  parts. 
By  using  the  symmetric  basis  weight  v(x)  =  [x(l  —  x)]3  and  a  full  complement  of  algebraic 
terms,  this  pitfall  can  be  avoided.  Furthermore,  the  spatial  basis  {(,}  as  defined  in  (3.36) 
can  be  used  for  problems  with  free  boundary  conditions,  thus  providing  consistency  to  the 
method. 

The  basis  for  the  problem  (3.2)  is  then  taken  to  be  {C«^}  with  Sj(t)  =  S(j,ht)  o  T(/), 
and  the  approximate  solution  is  defined  to  be 


«m.m,(*,0  =  E  £  «y<i(*)SJ(0 

i=-M§  j=-Mt 
Nt 

+  5i(0{u-W.-3,jC-M.-3(*)  +  U_Af._4,jC-M.-4(*) 

j=-Mt 

+UjV.+l,i(jV.+l(*)  +  tijv.+2,jCtf.+2(*)} 


+{a(0C-M.-i(*)  +  /?(0<-M.-2(*)  +  7(0C/v.+3(z)  +  £(t)(*.+4(s)} 

(3.38) 


where 


7(0  = 


1 

EI(1) 


7(0 


and 

=  Elflf®  ~  [£J(1)]*^‘ 

The  functions  7 (t)  and  6(t)  are  well-defined  since  EI(x)  is  assumed  positive  on  [0,1].  It 
should  be  noted  that  the  approximate  solution  does  satisfy  the  boundary  conditions  in  (3.2). 

The  (m*  +  4)  •  mt  unknowns  {u^}  in  (3.38)  are  determined  by  orthogonalizing  the  resid¬ 
ual  with  respect  to  the  sine  functions  {5p(x)5*(0}p=Im1-  v!,/v.+2-  This  Petrov-Galerkin 
approach  is  in  contrast  to  those  Galerkin  methods  in  which  the  residual  is  orthogonalized 
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with  respect  to  the  basis  and  is  done  to  take  advantage  of  the  exponential  accuracy  of  point 
evaluation  in  the  quadrature.  This  yields  the  discrete  system 


/,  SpSq)  —  0, 


-Mm  -2<p<Nm  +  2 
-Mt  <q<Nt 


where  (•,•)  is  defined  in  (3.5)  with  u»(x,t)  =  (T(t))-i. 

Appropriate  integration  by  parts  and  application  of  the  sine  quadrature  rules,  as  discussed 
in  [18],  yields  the  matrix  equation 


A(EI)UCj  +  CtMUAj  =  G 

where  CUCX,  and  At  are  defined  in  (3.20),  (3.21),  and  (3.23)  respectively,  and 


The  (mx  +  4)  x  mx  matrices  U  and  7  are  defined  componentwise  by 

[U)ij  =  No¬ 


where 


(7|.,  = 


7(X.  t)  s  /(*,  t)  -  £(3(i)B_„._,(*))  -  £(5(i)fl-M.-.(x)) 


-£(7(<)S«.«(x))  - 

The  (mx  -f  4)  x  (mx  +  4)  matrix  M  has  the  form 


M=  L 


where  the  mx  x  mx  submatrix  Dm  and  the  (mx  +  4)  x  1  vectors  are  given  by 


(3.39) 


Dm  —  D(v), 
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&L2  = 

3)1, 

&fll  =  'B{Blfm+ i)l 
and 

km  =  V(Bffm+ a)l. 

Here  1  i*  simply  an  (mx  +  4)  x  1  vector  of  ones. 

The  spatial  matrix  A(EI )  can  be  constructed  as  follows.  Let  and  p*(<)  be  defined 

as  they  were  in  (3.26),  (3.27),  and  (3.25),  respectively  (with  i  =  0, 1, 2  and  y  =  2,3, 4).  Note 
that  in  the  definitions  now,  the  index  ranges  are  — Mx  <  i  <  Nx,  —Mx  —  2<p<Nx  +  2,  and 
—Mx  —  2  <  k  <  Nx  +  2,  and  the  (mx+4)  x  1  coefficient  vector  is  now  c  =  [c_ ac«— a,  •  •  • ,  Cn.+a]T* 
Hence  and  p4(<)  have  the  sizes  (mx  +  4)  xmx,m.  x  (mx+4)  and  mx  x  1,  respectively. 

Furthermore,  let  4"  denote  the  (mx  +  4)  x  mx  matrix  which  is  defined  componentwise  by 

- wh(SrW™ 

and  let  c  =  {c_m..  •  •  • ,  cjy.]T.  Finally,  for  i  =  —  Mx— 4,  •  •  • ,  —  Mx— 1  and  i  =  Nx+l,  •  •  • ,  ATx+4, 
let  Si  denote  the  (mx  +  4)x  1  vectors 

*"V{veB';)c  +  V  ^(£/c£;r)  1]  • 

Here 

EIc{x)  =  C-Mu-\b-Mm- 1(*)  +  C_Af._j5_M._2(l)  +  Cffm+ibffm+l(x)  +  Cffm+2bfrm+l(x) 

and  1  is  simply  the  mx  x  1  vector  of  ones.  The  (mx  +  4)  x  (mx  -f  4)  matrix  A(EI )  is  then 

A(EI)  =  (o_m.-4  :  :  :  “W.+l  :  ®AT.+j]  (3.40) 

where  the  (mx  +  4)  x  mx  submatrix  Am  is  given  by 

Am  =  [$<a>P(u)Z>(p*(I))  +  2$<3>P(v)P(p*(i))  +  *<4>P(v)P(p*(o))] . 

It  should  be  noted  that  the  coefficient  matrix  A(EI)  in  (3.40)  differs  from  that  arising 
in  the  fixed  boundary  problem,  (3.24),  only  in  the  presence  of  the  diagonal  multipliers  V(v) 


21 


and  the  addition  of  border  vectors.  Hence  the  method  is  easily  adapted  when  changing  the 
boundary  conditions.  Furthermore,  the  matrices  4",  and  can  be  expanded  in  terms 
of  fundamental  matrices  in  a  manner  similar  to  that  in  (3.28)  -  (3.34),  thus  simplifying  the 
implementation  of  the  method. 

With  A(EI ),  At>  Cx,  Ct,  Af,  and  G  thus  specified,  the  system  (3.39)  can  be  solved  via  the 
generalized  Schur  algorithm  (3.35). 

The  final  implementation  issue  for  the  forward  problem  (3.2)  concerns  the  choice  of  decay 
parameters  a,/3, 7,  and  S.  As  discussed  in  [18],  the  weight  choice  w(x,t)  =  (T(t))“l  yields 
the  decay  condition 

\£I(x)U(x,t)\  <  *a+3(l  -  1  )*+V+ie-“  (3.41) 

where  £T{x)  is  defined  in  (3.11)  and  U(x,  t)  is  that  part  of  the  true  solution  which  is  ap¬ 
proximated  by 

«*(*.<)  =  £  13  «i,C<(*)5,-(<) 

(this  can  be  formally  obtained  by  subtracting  all  boundary  contributions  from  the  true 
solution  u(x,  <)).  With  the  decay  parameters  specified  and  Mm  chosen,  the  remaining  stepsizes 
and  summation  limits  are  given  by  (3.12)  -  (3.16). 


4  The  Finite  Dimensional  Minimization  Problem 


As  noted  in  the  introduction,  the  minimization  problem 


min  TJEI) 

EltQ  v  ' 


where 


1 


ra(Ei)=u\\ic(Ei)-dr+a\\EirQ} 


is  infinite  dimensional  and  thus  must  be  replaced  by  a  sequence  of  finite  dimensional  problems 
before  a  viable  numerical  scheme  can  be  developed.  Following  from  (3.10),  the  approximating 
admissible  parameter  sets  are  taken  to  be 


-l 


Nb 

Qmg  =  {  EIm„  :  EIma(x)  =  £  CfcV’fc(x)  J ,  mE  =  Me  +  Ne  +  1 

k=-Ma 


f 
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with  the  basis  {V>*}  defined  in  (3.8).  The  summation  limits  depend  on  the  boundary  condi¬ 
tions  with  Me  =  Mx,  Ng  =  Nx  for  fixed  boundary  conditions  and  Mg  =  Mx+2t  Ng  =  Nx+2 
for  cantilever  boundary  conditions.  The  associated  finite  dimensional  optimization  problem 
can  then  be  formulated  as 

<41> 

where 

=  i{||/t(E/„.)  -  If  +  a||J5/„X}-  (4.2) 

Note  that  in  solving  the  minimization  problem  (4.1),  one  is  actually  solving  for  the  vector 
c  =  [c_ •  •  • ,  catb]t  €  IR7'n  which  minimizes  Ta. 

With  rip  and  n,  specifying  the  number  of  spatial  and  temporal  observation  points,  re¬ 
spectively,  the  approximation  K(EImm)  :  MmM  -*  ]Rn*  n*  to  )C(EI)  is  obtained  by  applying 
a  discrete  analogue  of  the  observation  operator  C  in  (1.6)  to  Ummm,  in  (3.3)  or  (3.38).  If  the 
set  of  observation  points  {(*P,<9)}p=i|!!!|"j  can  be  represented  as  a  tensor  product  of  spatial 
and  temporal  points,  then  K(EIma)  has  the  representation 

k{EIm„)  =  C  co{U)  (4.3) 

where  the  matrix  U  solves  either  (3.19)  or  (3.39)  depending  on  the  boundary  conditions.  For 
the  fixed  boundary  problem  (3.1),  the  matrix  concatenation  co(U)  is  the  vector  in 
which  is  obtained  by  successively  stacking  the  columns  of  the  mx  x  mt  matrix  U.  C  is  an 
(rip-ti,)  x  (mx  mt)  evaluation  matrix  which  can  be  formulated  as  follows.  Define  the  np  x  mx 
spatial  evaluation  matrix  Ex  to  have  components 

[Ez]p,i  =  -Sifcp),  1  <  P  <  nP,  -Mx  <i<  Nx 

and  let  the  nq  x  mt  temporal  evaluation  matrix  Et  have  the  components 

[Et]q,j  =  l<q<nq,  -Mt<j<Nt. 


Then 


C  =  Ef  ®  Ex. 


Similar  formulations  for  C  and  co(U)  can  be  used  if  the  cantilever  boundary  value  problem 
is  being  considered.  It  is  noted  that  if  the  set  of  observation  points  is  not  rectangular  as 
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described  above,  then  point  evaluation  can  be  done  directly  via  (3.3)  or  (3.38).  This  latter 
option  is  less  efficient,  however,  than  that  defined  in  (4.3). 

The  discrete  penalty  functional  ||257m,|Jg  is  formed  by  substituting  the  expansion  (3.10) 
into  the  definition  (1.8).  This  yields 

||£/„.|lq  =  +  t  f[EIm,(x))'ix  a  lTQc 

JO  Jo 

where  the  mg  x  mg  matrix  Q  =  Qd  +  Qf  has  components 

[Qd]kt  «  [  ^k(x)^i{x)dx,  -Me  <  kyl  <  Ng 

JO 

and 

[Qf]ki  5=8  /  Mx)Mx)dx>  -ME<Kt<NE. 

JO 

The  matrix  entries  are  approximations  in  the  sense  that  sine  quadrature  rules  are  used  to 
evaluate  many  of  the  integrals. 

For  the  choice  of  basis  functions  in  (3.8),  the  matrix  Qd  is  given  by 

46  6  -62 

6  12  6  -12  6 

Qd=  0  0  Qd  6  6 

-6  -12  0  12  -6 

26  0  -64 

Integration  by  parts  and  the  application  of  the  sine  quadrature  formula  (2.11)  yields  the 
(mg  —  4)  x  (mg  —  4)  submatrix 


^  -  >  -  §> + r/°' 


where  again,  1^,1  =  0,2,4,  denote  the  matrices  whose  pi-th  entries  are  from  (2.12), 
(2.14),  and  (2.16),  respectively.  The  zeroing  of  all  other  quadrature  terms  is  a  result  of  the 
choice  vg( *)  =  (<£'(*))"$  =  [x(l  —  x)]i.  The  notation  0  simply  denotes  a  zero  vector  of 
length  mg  -  4.  Because  7(4)  is  positive  definite  and  7(2)  is  negative  definite  (see  [17]),  the 
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matrix  Q4  is  nonnegative  definite  with  zero  eigenvalues  resulting  from  the  configuration  of 
the  outer  four  columns. 

Direct  integration  and  sine  quadrature  are  used  to  obtain  the  matrix 


1 

11 

nT 

13 

-1 

105 

210 

9  12 

420 

140 

11 

13 

nT 

9 

13 

210 

35 

9  /l 

70 

420 

qt2 

911 

A 

Qf 

9r\ 

9t2 

13 

9 

r*T 

13 

11 

420 

70 

9  rl 

35 

120 

-1 

13 

=*T 

11 

1 

140 

420 

9  r2 

120 

105 

Here 

Qf  =  hxV(x\l  -  x)3) 

where  T>{tj)  denotes  the  diagonal  matrix  with  entries  tj(x_mb+2),'  •  •  Recall  that 

the  sine  gridpoints  are  defined  in  (2.7).  The  vectors  qti,qi2, qT\  and  qr2  have  components 

[£ilfc  =  Axx*(l  -  xfc)3[2xfc  -f  1J, 

\qa\k  =  hx(  1  -  xfc)xfc[2(l  -  xfc)  +  1], 

[?rl]fc  =  hxX2k{l  -  Xfc)3, 

and 

(9ra]fc  =  -hx(l  -  xfc)Vfc 

for  k  =  —Me  +  2,  •  •  • ,  Ne  —  2.  The  matrix  Qf  is  strictly  positive  definite. 

Although  the  matrix  Q  is  full,  it  is  very  efficient  to  construct  since  the  Toeplitz  matrices 
J(°),  and  1^  are  also  needed  in  the  forward  solver.  For  e  >  0,  Q  is  symmetric  and  positive 
definite  and  hence  has  a  Cholesky  decomposition  Q  =  RrR  where  R  is  upper  triangular.  It 
then  follows  that  the  penalty  term  ||£/fl,j||Q  yields  the  quadratic  form 

ctRtRc=  ||i?c  ||3  (4.4) 

where  j|  • )]  denotes  the  Euclidean  norm.  This  representation  for  the  penalty  functional  is 
particularly  useful  both  when  implementing  a  scheme  to  solve  the  minimization  problem  and 
when  plotting  the  L-curve  to  determine  a  suitable  regularization  parameter  a  (see  [6]). 
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To  highlight  the  dependence  of  the  functional  Ta  in  (4.2)  on  the  unknown  vector 


c  =  [c-mb>  •  •  • ,  cjvs]t,  let 


K(c)  =  )  =  C  co{U  (c )) 


where  U(c )  solves  either  (3.19)  or  (3.39)  for  a  given  expansion  EIma(x)  = 
Noting  (4.4),  the  optimization  problem  (4.1)  can  be  replaced  by 


nb 

£  ckxf)k(x). 
k-~Mg 


min 

^€Rm* 


US) 


(4.5) 


where 

US)  =  jdl*  (?)  -  <*1l’  +  «»«?«’}•  (4-6) 

To  obtain  a  minimizer  for  the  nonlinear  functional  Ta(c),  a  quasi-Newton/trust  region  iter¬ 
ation  [2] 

Cfc+i  =  ck  +  sk 


is  used.  Here  sk  solves  the  quadratic  programming  problem 

,min.  \  {ll«-(5>)  +  K'(S„)S~  if  +  aHm  +  ?)l|!} 

subject  to  ||s  ||  <  6k  with  K'{ck )  denoting  the  Jacobian  of  K  at  ck.  The  trust  region  radius 
8k  is  chosen  so  that  Ta(c )  has  sufficient  decrease  at  each  iteration  to  guarantee  convergence 
to  a  local  minimizer  of  Ta  (for  further  details  about  the  theory  and  implementation  of  the 
trust  region  algorithm,  see  [2]  or  [8]).  An  important  numerical  issue  in  the  implementation 
of  the  trust  region  scheme  is  the  formulation  of  the  derivative  of  the  operator  K.  Here  the 
derivative,  or  Jacobian,  is  an  (np  •  n9)  X  matrix  whose  v-th  column  is  given  by 

\K\S) ]„  =  limi|K(c'+  TK)  -  *(c)l 
where  the  standard  unit  vector  e„  has  components 


[2|/]fe  =  f>vk  =  { 


1 


if  k  =  «/, 


—Me  <  k  <  Ne 


y  0  otherwise. 

In  the  examples  of  the  next  section,  the  Jacobians  were  calculated  with  a  standard  forward 
difference  scheme.  This  scheme  is  easy  to  implement  and  accurate  enough  for  the  purposes 
of  the  method.  If  further  efficiency  is  desired,  a  directional  derivative  scheme  such  as  that 
described  in  [12]  can  be  used. 
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5  Implementation  and  Numerical  Examples 


The  four  examples  reported  in  this  section  were  selected  from  a  large  collection  of  problems 
to  which  the  Sinc-Galerkin  method  was  applied.  The  results  are  representative  of  those 
obtained  for  other  problems. 

The  first  example  demonstrates  the  application  of  the  Sinc-Galerkin  method  to  a  model 
problem  with  fixed  boundary  conditions  in  which  the  state  solution  was  sampled  directly; 
that  is,  no  external  noise  was  added  to  the  data.  To  demonstrate  the  feasibility  of  the  method 
for  problems  with  noisy  data,  the  same  problem  is  revisited  in  Example  5.2  but  with  pseudo¬ 
random  white  noise  added  to  the  data.  In  the  third  example,  the  parameter  to  be  recovered 
is  the  shifted  Gaussian  function  that  was  discussed  for  second-order  examples  in  [12]  and  [19]. 
Again,  pseudo-random  white  noise  is  added  to  the  data.  The  final  example  demonstrates 
the  application  of  the  method  to  a  problem  with  cantilever  boundary  conditions  as  modeled 
by  (1.2).  Hence  in  this  example,  the  parameter  El  appears  both  in  the  Bpatial  operator  and 
in  the  boundary  conditions. 

In  all  examples,  d  =  |  (see  (2.4)  and  (2.6)).  The  errors  for  the  method  are  reported  on 
the  set  of  uniform  gridpoints 

U  =  {zo,Zi,  ’ '  izioo) 
zk  =  fcl,  k  =  0, 1,-  •  • ,  100 

with  stepsize  t  =  -—g.  With  El  and  EIma  denoting  the  true  and  approximate  parameters 
respectively,  the  errors  are  reported  as 

IIM<)II  =  (5.2) 

The  error  results  are  tabulated  in  the  form  .ooo  —  7  which  represents  .aaa  x  lO-'*  and  all 
problems  were  run  with  sixteen  place  accuracy  on  a  Vax  8550. 

A  very  important  practical  consideration  is  the  choice  of  the  regularization  parameter  a 
for  a  given  (error-contaminated)  data  set.  If  the  error  in  the  data  is  discrete  and  random, 
then  under  certain  conditions  the  method  of  Generalized  Cross  Validation  (GCV)  can  be  used 
to  determine  a  suitable  value  of  a  [25].  A  second  method  for  determining  the  regularization 
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parameter  is  to  plot  the  norm  of  the  penalty  functional,  ||Aca||,  versus  the  norm  of  the 
residual,  ||A(ca)  —  <f||  (see  [6]  or  [14]).  Here  ca  denotes  the  solution  to  (4.5).  In  this  way, 
one  can  qualitatively  get  an  idea  of  the  compromise  between  the  minimization  of  these  two 
quantities.  The  scheme  for  determining  the  “optimal”  regularization  parameter  consists  of 
finding  those  values  of  a  such  that  (||A(ca)  —  <f||,  ||i2ca||)  lies  in  the  “corner”  of  the  resulting 
curve,  known  as  the  X-curve.  The  use  of  this  technique  for  determining  suitable  choices  for 
the  regularization  parameter  is  demonstrated  in  the  examples. 

In  all  four  examples,  the  data  was  sampled  on  a  regular  grid  {(xp,  $,)}  C  (0, 1)  x  (0,2]. 
Nineteen  equally  spaced  points  xp  =  pAx,  Ax  =  .05,  were  taken  in  space  and  four  equally 
spaced  temporal  points  tq  =  qAt,  At  =  .5,  were  taken  for  a  total  of  n  =  76  data  points.  In 
all  examples,  the  mg  x  1  initial  vector  3)  =  [.5,  .5,  .5,  •  •  • ,  .5,  .5,  .5]r  was  used. 

Finally,  it  should  be  noted  that  in  the  examples,  the  symbol  a  is  used  to  denote  both 
the  regularization  parameter  and  the  sine  decay  parameter.  The  use  of  this  symbol  for  both 
quantities  is  well  established  in  the  literature  and  is  thus  difficult  to  avoid  in  this  setting. 
It  should  be  obvious  from  the  context  however,  which  quantity  is  being  discussed  and  there 
should  be  no  ambiguity  resulting  from  the  dual  use  of  this  symbol. 


Example  5.1. 


d7u 
Tt T 


d 2  /  d2u\ 

+  w[EIWw)  =  «*■*■  t>0  0 

u(0,  t)  =  tt(l,  t)  =  0,  t>0 

!«.,<)= ^(1,0=0,  «>o 

du 

«(*, o)  =  ^-(x,0)  =  0  <  X  <  1 


<  x  <  1 


The  forcing  function  /(x,  t)  is  consistent  with  the  true  stiffness  parameter  EI(x)  =  l+sin(irx) 
and  the  state  solution  u(x,t)  =  x(l  —  x)sin(4irx)tae-‘.  For  these  functions,  the  choices 
a  =  (i  =  7  ss  |  and  6  =  1  satisfy  the  decay  condition  (3.17).  No  noise  was  added  so 
the  data  consisted  of  direct  measurements  of  the  state  solution.  For  varying  values  of  the 
regularization  parameter  a,  the  L- curve  is  plotted  in  Figure  3.  Note  that  the  value  a  =  10“8 
yields  a  point  (||A(ca)  -  <f||,  ||i?ca||)  in  the  "corner”  of  the  curve.  The  uniform  errors  for 
a  =  10“8  are  reported  in  Table  1  with  the  first  four  columns  indicating  the  index  limits 
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for  the  expansion  of  the  state  variable  and  the  fifth  column  indicating  the  number  of  basis 
functions  used  in  the  expansion  of  EIma.  The  convergence  of  the  method  is  demonstrated 
both  by  the  results  in  the  last  column  of  Table  1  and  by  Figure  4  which  shows  the  true  and 
approximate  stiffness  parameters  with  a  —  10-8. 


Ns 

Mt 

Nt 

mE 

\\EIu(l)\\ 

8 

8 

8 

4 

17 

.1869  -  0 

16 

16 

16 

6 

33 

.2482  - 1 

24 

24 

24 

7 

49 

.1463  -  2 

Table  1.  Errors  on  the  Uniform  Grid  U  with  a  =  10“8  in  Example  5.1. 


Figure  3.  The  Tikhonov  Z-Curve  for  Example  5.1. 


Figure  4.  True  and  Approximate  Stiffness  Parameters  for  Example  5.1  with  a  =  10~B 
-  •  -  (mE  =  17), - (mE  =  33),  •  •  •  (mg  =  49),  - (True). 

Example  5.2. 

In  this  example,  the  true  parameter  and  state  solution  are  the  same  as  those  in  Example  5.1, 
and  hence  EI(x)  =  1  +  sin(Trx)  and  u(x,t)  =  x(l  -  x)sin(4xx)t3e_*.  To  the  data  however, 
we  added  a  pseudo-random  noise  vector  e  from  a  Gaussian  distribution  with  mean  0  and 
standard  deviation  <r  chosen  so  that  the  noise- to-signal  ratio  o-/|(d  ||  =  .01  ;  that  is,  noise 
=  1%  of  the  signal.  The  L-curve  is  plotted  in  Figure  5.  Note  that  the  values  a  =  10-B  and 
a  =  5  x  10-8  yield  points  (||/f(c^)  —  d  ||,  ||/2ca||)  in  the  “corner"  of  the  curve.  For  mg  =  33, 
the  uniform  errors  obtained  with  a  =  10“3,  a  =  10-5  and  a  =  10-1°  are  given  in  Table  2. 
Corresponding  plots  of  the  true  and  approximate  stiffness  parameters  are  shown  in  Figure  6. 
Note  that  the  “corner”  value  a  =  10-8  provides  a  very  good  choice  for  the  regularization 
parameter  whereas  a  =  10_1°  is  not  large  enough  to  damp  out  the  error  contributions  due 
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to  the  smaller  singular  values.  Finally,  the  choice  a  =  10-3  causes  too  much  smoothing  and 
information  about  the  parameter  is  lost.  The  results  from  this  example  demonstrate  the 
viability  of  the  method  for  problems  with  noisy  data. 


a  =  10~3  a  =  10“8  a  =  10“10 

\\EIu(t)\\ 

.8503  -  0  .2228  -  1  .3840  -  1 

Table  2.  Errors  on  the  Uniform  Grid  U  with  mj  =  33  in  Example  5.1. 


Residual  Norm  |Jf(&)  -  / R 

Figure  5.  The  Tikhonov  L-Curve  for  Example  5.2. 
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2.2 


Ofli  • _ i  i _ i _ i _ i _ i - ■ - ■ - 1 

0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  I 

x-ixi* 


Figure  6.  True  and  Approximate  Stiffness  Parameters  for  Example  5.2  with  mg  =  33 
- (a  =  10-3), - (a  =  10“8),  •••(«=  10~10),  - (True). 


Example  5.S. 


W+&{EI{x)l&)=K*’t)’  <>0  0<x<1 

u(0,<)  =  u(l,i)  =  0,  t  >  0 

§><)=fj(M)  =  o,  <>o 

«(*.  0)  =  ^(*,0)  =  0,  0  <  x  <  1 

In  this  example  the  stiffness  parameter  to  be  recovered  is  the  shifted  Gaussian  function 
EI(x)  —  1  +  The  state  solution  u(x,t)  =  sin3(xx)iJe“‘  yields  the  decay  param¬ 

eters  a  =  0  =  7  =  |  and  6  =  1  as  dictated  by  (3.17).  Pseudo-random  noise  is  again  added 
to  the  data  in  the  manner  described  in  the  last  example.  As  seen  in  Figure  7,  the  Tikhonov 
parameter  values  a  =  10“*  through  a  =  5  x  10“8  yield  points  (||/C(5x)  -  A  ||,  ||J2c*||) 
in  the  “corner”  of  the  L-curve.  For  the  “corner”  value  o  =  10“7,  numerical  results  with 
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mg  =  17,  ms  =  33  and  mg  =  49  are  reported  in  Table  3  and  Figure  8.  In  spite  of  the  noise 
in  the  data,  both  sets  of  results  demonstrate  that  the  method  converges  as  the  number  of 
basis  functions  is  increased. 


Af, 

Nt 

Mt 

Nt 

mg 

\\EIu{t)  || 

8 

8 

8 

4 

17 

.1811-0 

16 

16 

16 

6 

33 

.4191  -  1 

24 

24 

24 

7 

49 

.1058  - 1 

Table  3.  Errors  on  the  Uniform  Grid  U  with  a  =  10  7  in  Example  5.3. 


Figure  7.  The  Tikhonov  L-Curve  for  Example  5.3. 
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Figure  8.  True  and  Approximate  Stiffness  Parameters  for  Example  5.3  with  a  =  10  7 
- (mE  =  33), - (mE  =  49),  - (True). 

Example  5-4 . 

t>0  0<x<1 

u(o,()  =  o,  (i,0  =  o.  t  >  o 

|(o.o  =  1V<,  £(*/£)< l.D-faVr*.  oo 

a 

u(i,0)=^(*,0)  =  0,  0<I<1. 

This  example  demonstrates  the  Sinc-Galerkin  method  for  a  problem  with  cantilever  boundary 
conditions  and  hence  the  stiffness  parameter  appears  both  in  the  spatial  operator  and  in  the 
boundary  conditions  themselves.  The  true  stiffness  parameter  is  EI(x)  —  l+sin(xx)  and  the 
state  solution  is  u(z,<)  =  sin3(x®)tJe“‘.  For  these  functions,  the  choices  a  =  /?=l,7=§ 
and  S  =  1  satisfy  the  decay  condition  (3.41).  No  noise  was  added  so  the  data  consisted 
of  direct  measurements  of  the  state  solution.  Since  the  L-curve  was  very  similar  to  that  of 
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Example  5.1,  the  regularization  parameter  was  taken  to  be  a  =  10~8.  The  uniform  errors 
for  this  choice  are  reported  in  Table  4  and  the  true  and  approximate  parameters  are  shown 
in  Figure  9.  When  comparing  Tables  4  and  1,  it  is  noted  that  the  index  limits  differ  as  a 
result  of  the  choice  a  =  /3  =  1  for  the  decay  parameters.  The  smaller  values  of  a  and  /3  also 
indicate  why  the  errors  here  are  slightly  larger  than  those  in  Example  5.1.  Finally,  because 
Mg  =  Mx+2  and  Ne  =  Nx  +  2  for  these  boundary  conditions,  the  number  of  basis  functions, 
ms,  used  in  the  expansion  of  EIma  also  differs  from  the  number  used  in  Example  5.1  where 
Me  =  Mx  and  Ne  =  Nx.  Both  the  table  and  the  figure  demonstrate  the  convergence  of  the 
method  for  problems  with  cantilever  boundary  conditions. 


Ms  N,  Mt  Nt  mE  \\EIv(i)\\ 


8 

8 

6 

3 

21 

.1589 

-0 

16 

16 

11 

4 

37 

.2589 

-1 

24 

24 

16 

6 

53 

.1334 

-  1 

Table  4.  Errors  on  the  Uniform  Grid  U  with  a  =  10~8  in  Example  5.4. 


Figure  9.  True  and  Approximate  Stiffness  Parameters  for  Example  5.4  with  a  =  10 
-  •  -  (ms  =  21), - (ms  =  37),  •  •  •  (mg  =  53),  - (True). 


References 


[1]  E.A.  Coddington  and  N.  Levinson,  Theory  of  Ordinary  Differential  Equations , 
McGraw-Hill,  New  York,  1955. 

[2]  J.E.  Dennis  and  R.B.  Schnabel,  Numerical  Methods  for  Unconstrained  Optimization 
and  Nonlinear  Equations ,  Prentice-Hall,  Englewood  Cliffs,  NJ,  1983. 

[3]  G.H.  Golub,  S.  Nash  and  C.  VanLoan,  “A  Hessenberg-Schur  Method  for  the  Prob¬ 
lem  AX  +  XB  =  C,”  IEEE  Trans.  Automat.  Control,  24,  909-913  (1979). 

[4]  G.H.  Golub  and  C.  VanLoan,  Matrix  Computations,  2nd  ed.,  Johns  Hopkins  Uni¬ 
versity  Press,  Baltimore,  1989. 

[5]  R.  Guenther  and  J.  Lee,  Partial  Differential  Equations  of  Mathematical  Physics 
and  Integral  Equations ,  Prentice-Hall,  Englewood  Cliffs,  NJ,  1988. 

[6]  P.C.  Hansen,  “Analysis  of  Discrete  Ill-Posed  Problems  by  Means  of  the  L- curve,” 
submitted  to  SIAM  Rev. 

[7]  K.  Jonca,  "Numerical  Solution  of  a  Nonlinear  Fredholm  Integral  Equation  of  the 
First  Kind,”  PhD  Dissertation,  Montana  State  University,  1988. 

[8]  K.  Jonca  and  C.R.  Vogel,  “Numerical  Solution  to  the  Magnetic  Relief  Problem,” 
in  Transport  Theory,  Invariant  Imbedding,  and  Integral  Equations  (Lecture  Notes 
in  Pure  and  Applied  Mathematics )  ed  P  Nelson  et  al,  Marcel  Dekker,  New  York, 
379-391  (1989). 

[9]  K.  Kunisch  and  E.  Graif,  “Parameter  Estimation  for  the  Euler- Bernoulli  Beam,” 
Mat.  Aplicada  e  Computational,  4,  95-124  (1985). 

[10]  W.  Leighton  and  Z.  Nehari,  “On  the  Oscillation  of  Solutions  of  Self-Adjoint  Linear 
Differential  Equations  of  the  Fourth  Order,”  Trans.  Am.  Math.  Soc.,  89,  325-377 
(1968). 

[11]  J.  Lund,  “Symmetrization  of  the  Sinc-Galerkin  Method  for  Boundary  Value  Prob¬ 
lems,”  Math.  Comp.,  47,  571-588  (1986). 


[12]  J.  Lund  and  C.R.  Vogel,  “A  Fully-Galerkin  Method  for  the  Solution  of  Inverse 
Problems  in  Parabolic  Partial  Differential  Equations,”  Inv.  Prob.,  6,  205-217  (1990). 

[13]  K.M.  McArthur,  K.L.  Bowers  and  J.  Lund,  “Numerical  Implementation  of  the 
Sinc-Galerkin  Method  for  Second-Order  Hyperbolic  Equations,”  Numer.  Methods 
Partial  Diff.  Equations,  3,  169-185  (1987). 

[14]  K.  Miller,  “Least  Squares  Methods  for  Ill-Posed  Problems  with  a  Prescribed 
Bound,”  SIAM  J.  Math.  Anal.,  1,  52-74  (1970). 

[15]  C.B.  Moler  and  G.W.  Stewart,  “An  Algorithm  for  Generalized  Matrix  Eigenvalue 
Problems,”  SIAM  J.  Numer.  Anal.,  10,  241-256  (1973). 

[16]  T.I.  Seidman  and  C.R.  Vogel,  “Well-Posedness  and  Convergence  of  Some  Regular¬ 
ization  Methods  for  Nonlinear  Ill-Posed  Problems,”  Inv.  Prob.,  5,  227-238  (1989). 

[17]  R.C.  Smith,  G.A.  Bogar,  K.L.  Bowers  and  J.  Lund,  “The  Sinc-Galerkin  Method 
for  Fourth-Order  Differential  Equations,”  to  appear  in  SIAM  J.  Numer.  Anal. 

[18]  R.C.  Smith,  K.L.  Bowers  and  J.  Lund,  “A  Fully  Sinc-Galerkin  Method  for  Euler- 
Bernoulli  Beam  Models,”  submitted  to  Numer.  Methods  Partial  Diff.  Equations. 

[19]  R.C.  Smith  and  K.L.  Bowers,  “A  Fully  Sinc-Galerkin  Method  for  Parameter  Esti¬ 
mation  in  Parabolic  Problems,”  to  be  submitted  to  Inv.  Prob. 

[20]  I.  Stakgold,  Green’s  Functions  and  Boundary  Value  Problems,  Wiley,  New  York, 
1979. 

[21]  F.  Stenger,  “Approximations  Via  Whittaker’s  Cardinal  Function,”  J.  Approx.  The¬ 
ory,  17,  222-240  (1976). 

[22]  F.  Stenger,  “A  Sinc-Galerkin  Method  of  Solution  of  Boundary  Value  Problems,” 
Math.  Comp.,  33,  85-109  (1979). 

[23]  F.  Stenger,  “Numerical  Methods  Based  on  Whittaker  Cardinal,  or  Sine  Functions,” 
SIAM  Rev.,  23,  165-224  (1981). 


37 


(24]  A.N.  Tikhonov  and  V.Y.  Arsenin,  Solutions  of  Ill-Posed  Problems ,  Wiley,  New 
York,  1977. 

[25]  G.  Wahba,  “Practical  Approximate  Solutions  to  Linear  Operator  Equations  when 
the  Data  are  Noisy,”  SIAM  J.  Numer.  Anal.,  14,  651-667  (1977). 


38 


